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A health monitoring system based on analytical redundancy is developed for satellites 
on elliptical orbits. First, the dynamics of the satellite including orbital mechanics and at- 
titude dynamics is modelled as a periodic system. Then, periodic fault detection filters are 
designed to detect and identify the satellite’s actuator and sensor faults. In addition, parity 
equations are constructed using the algebraic redundant relationship among the actuators 
and sensors. Furthermore, a residual processor is designed to generate the probability of 
each of the actuator and sensor faults by using a sequential probability test. Finally, the 
health monitoring system, consisting of periodic fault detection filters, parity equations 
and residual processor, is evaluated in the simulation in the presence of disturbances and 
uncertainty. 


I. Introduction 

Satellites under automatic control demand a high degree of reliability in order to operate properly. If a 
sensor fails, the controller’s command will be generated using incorrect measurements. If an actuator fails, 
the controller’s command will not be applied properly to the satellite. To avoid this situation, one needs 
a health monitoring system capable of detecting a fault as it occurs and identifying the faulty component. 
This process is called fault detection and identification (or isolation). 

The most common approach to fault detection and identification is hardware redundancy which is the 
direct comparison of the outputs from identical components. It requires very little computation. However, 
hardware redundancy is expensive and limited by space and weight. An alternative is analytical redundancy 
which uses the modelled dynamic relationship between control commands and measurements to form a 
residual process. Nominally, the residual is nonzero only when a fault has occurred and is zero at other 
times. Furthermore, the residual is nonzero in a unique and fixed direction in response to each fault. 
Therefore, when a nonzero residual is detected, a fault can be announced and the faulty component can be 
identified. 

One popular approach to analytical redundancy is the restricted diagonal detection filter 1,2 which has 
two special cases. When every fault is detected, it becomes the Beard-Jones detection filter 3,4 which can 
be determined by spectral methods. 2, 5,6 When only one fault is detected, it becomes the unknown input 
observer. '~ 9 Ref. 10 to 12 generalize the unknown input observer for time-invariant systems to the approx- 
imate unknown input observer for time-varying systems. Furthermore, Ref. 13 generalizes the approximate 
unknown input observer in Ref. 12 from detecting single fault to multiple faults to obtain the approximate 
restricted diagonal detection filter for time- varying systems. 

In Ref. 14, the fault detection filters in Ref. 12 and 13 are applied to periodic systems. For the approximate 
unknown input observer, the filter gain becomes periodic naturally. For the approximate restricted diagonal 
detection filter, the filter gain becomes periodic by imposing a boundary constraint on the initial and final 
conditions. A numerical algorithm is given to solve for the periodic filter gain. This is an important extension 

* Research Engineer, Mechanical and Aerospace Engineering Department 
^ Graduate Student Researcher, Mechanical and Aerospace Engineering Department 

Professor, Mechanical and Aerospace Engineering Department, AIAA Fellow 
§ Graduate Student Researcher, Mechanical and Aerospace Engineering Department 
^Aerospace Engineer, Flight Dynamics Analysis Branch, Code 595, AIAA Senior Member. 


1 of 22 


American Institute of Aeronautics and Astronautics 



for the approximate restricted diagonal detection filter because the periodic filter gain can be computed off- 
line and stored to be used repetitively. 

Although the residual generated by the fault detection filter is nominally zero in the absence of a fault 
and nonzero otherwise, when driven by disturbances and uncertainty, the residual fails to go to zero even 
in the absence of a fault. Therefore, it is difficult to determine the occurrence of a fault by simply choosing 
a threshold for the residual. To enhance fault detection and identification, a residual processor is used to 
analyze the residual which can be viewed as a pattern containing information about the presence or absence 
of a fault. A hypothesis is determined according to the pattern which is unique in response to a particular 
fault. By considering the residual processor design as a hypothesis testing problem, the residual processor 
could be a multiple-hypothesis Shiryayev sequential probability test. 15 The fault identification problem is 
now solved by assuming that each fault corresponds to certain hypothesis. Therefore, the residual processor 
generates the conditional probability of each fault hypothesis rather than a binary announcement. This 
allows for higher level decision making which now can be based on the probability. 

In this paper, a health monitoring system based on analytical redundancy is developed for satellites 
on elliptical orbits. Different from Ref. 14 where only orbital mechanics is included and the orientation of 
the actuators and sensors is assumed to be fixed in the ECI reference frame, both orbital mechanics and 
attitude dynamics are now included and the orientation of the actuators and sensors is assumed to be fixed 
in the satellite reference frame to reflect a more realistic situation. Furthermore, in addition to periodic 
fault detection filters, the residual processor is used to develop a complete health monitoring system for 
satellites. In Section II, the background of the periodic fault detection filters and residual processor is given. 
In Section III, the satellite on an elliptical orbit is modelled as a periodic system. In Section IV, periodic fault 
detection filters are designed and evaluated in the presence of disturbances and uncertainty. In Section V, 
the residual processor is designed and evaluated using the residuals generated in Section IV. 


II. Background 

In Section A, the actuator and sensor fault model is given. In Section B, the periodic fault detection 
filters in Ref. 14 are given. In Section C, the multiple-hypothesis Shiryayev sequential probability test in 
Ref. 15 is given. 


A. Actuator and Sensor Fault Model 

Consider a linear periodic system, 

x = Ax + Bu (la) 

y = Cx (lb) 

where u is the control command, y is the measurement, A(t + T) = A(t), B{t + T) = B{t), C(t + T) = C(t) 
and T is the period. The ith actuator fault can be modelled as an additive term in the state equation (la) 

x = Ax + Bu + F a y, a 

where the fault direction F a is the ith column of B and the fault magnitude y a is an unknown and arbitrary 
scalar function of time that is zero when there is no fault. 

The ith sensor fault can be modelled as an additive term in the measurement equation (lb). 

y = Cx + E s n s (2) 

where the fault direction E s is a column of zeros except a one in the ith position and the fault magnitude 
y s is an unknown and arbitrary scalar function of time that is zero when there is no fault. For the purpose 
of fault detection filter design, an input to the state equation that drives the measurement in the same way 
that /j, s does in (2) is obtained. 10 Define a new state x = x + where E s = Cf s . Then, the dynamic 
equation of x and (2) can be written as 



(3a) 

(3b) 
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where f s = f s — Af s . Therefore, for fault detection filter design, the sensor fault is modelled as a two- 
dimensional additive term in the state equation as in (3). 

B. Periodic Fault Detection Filter 

The periodic multiple-fault detection filter problem considers a linear periodic system with q faults. 

Q 

x = Ax + B u u + B w w + FiHi 

i= 1 

■y = Cx + v 

where w is the process noise and v is the sensor noise. If only s of the q faults need to be detected and 
identified, the objective of the periodic multiple-fault detection filter problem is to find a periodic filter gain 
L for the linear filter, 


x = Ax + B u u + L(y — Cx) 
r = y — Cx 


such that each projected residual H t r for i = 1 • • • s is affected primarily by its associated target fault Hi and 
minimally by its associated nuisance fault fi t = [/./-] • • • m~l lM+i ■ • • /i 9 ] T , process noise w, sensor noise v 
and initial condition error x(to) — x(to)- The projector Hi is defined a priori 10 


Hi = I- Ker H, 


(Ker H t f Ker H, 


{KerHf 


where Ker H = [Cbi t i^g il Cbi t 2 ,g i2 ••• Cbi >Pi> g ip . ]■ The vectors bij^ i:j , j = 1 •••pu are found from the 
iteration defined by the Goh transformation, 16, 17 

bi.j h k = Abijk—i bij : k — l 
bi.jji fi,j 

where f. t / is the jth column of F, = [ Fj ■ ■ ■ F, l+t ■ ■ ■ F q ] and p, = dimFj. S. t ] is the smallest non- 
negative integer such that Cbij ! g i . 0. Note that the derivation of Hi can also be included in the periodic 

multiple-fault detection filter problem instead of being defined a priori. 

Here the optimization problem for the periodic multiple-fault detection filter is presented after some 
manipulations. Please refer to Ref. 14 and 13 for more details. 


min 

L,Wi(t 0 )-W s (t 0 ) 



^ HiCWiC T Hi 

_i = 1 


dt 


(4) 


subject to 

Wi = (A- LC)Wi + Wi(A - LCf + (L — PiC T V- l )V{L - P, ; C' 7 V- 1 ) T (5a) 

0 = W i (t o + T)-W i (to) (5b) 

where 

P, = APi + P t A T - PiC T V- x CPi + -FiQiFj - F.Q.Fj + B W Q W B T W ( 6 ) 

7 * 

for i = 1 • • • s. For each projected residual H,r, Qi and 'y~ 1 Qi are the design weightings on the transmission 
from the associated target fault and associated nuisance fault, respectively. When Qi is larger, the trans- 
mission from /Xi is larger. When 'y^ 1 Qi is larger, the transmission from /j, is smaller. Similarly, Q w , V and 
Pq are the design weightings on the transmission from the process noise, sensor noise and initial condition 
error, respectively. When Q w , V and Pq are larger, the transmission from w, v and x(to) — x(to) is smaller, 
respectively. Note that the Riccati matrix P, in (6) is periodic and independent of L. 
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By taking the first-order variation of (4) and (5), and using continuously differentiable matrix Lagrange 
multipliers Ki s to adjoin the constraint to the cost, the first-order necessary condition implies that the 
optimal solution for L and the dynamics of are 


L* 


Z K ‘ 


Z K i(Pi + Wi) 


cFv- 1 


( 7 ) 


1 


Li=l 


satisfying (5) and 


-Ki = Ki(A - LC) + {A - LC) T Ki + C T HiC 

0 = Ki(t 0 + T) - Ki(t 0 ) 


(8a) 

(8b) 


where i = 1 • • • s. The process of deriving the filter gain of the periodic multiple-fault detection filter requires 
the solution to a two-point boundary value problem that includes a set of Lyapunov equations, (5) and (8), 
coupled by (7). A numerical algorithm is given in Ref. 14 to solve (4) numerically by using the steepest 
descent method. However, the global minimum cannot be guaranteed because (4) may not be convex. 

For the periodic single-fault detection filter problem that detects the fault /q and blocks the fault the 
minimization problem (4) reduces to 


1 

mm 

L t\ — < 0 



HiCWiC T Hi 


cLt 


subject to (5). The optimal filter gain is 

L* = PiC T V~ l (9) 

using the solution to (6). Note that s single-fault detection filters are needed in order to monitor all s faults. 
For all the results in the application in Section IV, the periodic single-fault detection filter is used. 


C. Residual Processor 


The residual processing problem is considered as a hypothesis testing problem which can be solved by the 
multiple- hypothesis Shiryayev sequential probability test (MHSSPT). 15 MHSSPT is a generalized result of 
Ref. 18 which solves the change detection problem based on the results of Ref. 19 and by using a dynamic 
programming formulation. Under certain condition, MHSSPT detects the change in hypothesis (i.e. , the 
occurrence of a disruption) in a conditionally independent measurement sequence in minimum time. In the 
dynamic programming formulation, the measurement cost, the cost of false alarm and the cost of miss alarm 
are considered. MHSSPT is shown to be optimal in the infinite time case. 

Here, the propagation equation for the posterior probability of disruption conditioned on the measurement 
sequence will be introduced first. This recursive formula allows the posterior probability of each hypothesis 
to be calculated online and plays a central role in the dynamic programming analysis. Then, the assumptions 
behind the derivation are discussed when the recursive formula is applied to our application. The propagation 
equation of the posterior probability of each hypothesis based on the measurement history is 


where 


Fk+l,i 
P(0i < tk+ i/X k ) 

P{0o < tk+i/Xk) 

Fo,i 


P(9j < tk+i/X k ) ■ fi{x k + 1 ) 

EZo < tk+1 / X k ) ■ fiixk+i) 

Fk,i + Pi(l — Fk,i) Vi ^ 0 

ad - P & ^ tk+l/X k )) Vi yf 0 

7T,: 


(10) 

( 11 ) 

(12) 

(13) 


Fk,i '■ P{0i < t k /X k ) 

9i : Time of transition to hypothesis Ti, 

Xk : Measurement history up to tk 

pi : A priori probability of transition to hypothesis 7 V; from tk to ifc+i, Vfc 
fi(xk) ■ Probability density function of Xk given Tti 
m + 1 : Number of hypotheses 
7b : P(0i < t 0 ). 
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The propagation equation (10) is derived from Bayes Rule under the following assumptions. First, 
the measurement sequence Xk is conditionally independent or equivalently, P(Xk/0i < tk) = P{xk/0i < 
tk)P(xk~i/0i < tk) ■ ■ ■ P(xi/6i < tk) which means that the measurement sequence is assumed to be 
independent when a disruption occurred. Second, the statistical properties of the measurements (i.e., the 
probability density function /)(•)) are assumed known before and after disruption for all hypotheses. Third, 
the a priori transition probability pi is assumed to be known and constant for all stages. Note that the analysis 
remains the same even if the transition probability is stage dependent. Finally, the a priori probability 7 q is 
also assumed to be known for all hypotheses. 

Now the validity of the above assumptions are discussed when the recursive formula is applied to our 
application. First, it is likely that each of the residual process generated by the periodic fault detection 
filter is time correlated. In our application, the collection of residuals generated at each time instant is 
considered independent static pattern which can be associated to one of the pre-defined hypothesis. Second, 
the distribution of the measurement sequence is assumed to be known before and after disruption. In practice, 
the magnitude of the failure is not known priori. The statistical properties of the distribution defined for all 
failure hypotheses may not be completely known. To deal with this uncertainty, if one of the parameters a 
of the density function _/)(•) is not known and assumed to follow its own distribution (i.e., a density function 
ip a (') defined over fi), then the conditional density function /,(•) can be formulated as 

/»(•) = [ fi{x/ri) ■ ipaiv) ' dr] (14) 

J n 

Finally, the probability of transition pi and the a priori probability 7 q are not known. Both of them are treated 
as design parameters. In general, a smaller pi is chosen when noisy measurement sequence is processed. The 
a priori probability 7r,; is assumed to be a small number. These two parameters are not sensitive in the 
algorithm and should not affect the performance when reasonable values are assumed. 

III. Satellite Model 

In this section, an Earth pointing satellite, which keeps one face of the satellite pointed to the Earth 
all the time, is modelled as a periodic system. The actuators considered are three thrusters. The sensors 
considered are a star sensor, a horizon sensor, a sun sensor, a GPS and an inertial measurement unit which 
has three gyros and three accelerometers. In Section A, the dynamics of the satellite is given. In Section B, 
the measurement model is given. In Section C, the dynamics of the satellite and the measurement model 
are linearized for the periodic fault detection filter design. 

A. Dynamics of Satellite 

The equation of motion associated with the orbital mechanics of the satellite is 

x + ||^j| 3 r — CgC^a = 0 (15) 

where r = [r x r y r z ] T is the position of the satellite in ECI reference frame, a = [a x a y a z } T is the accel- 
eration generated by the three thrusters in the thruster reference frame, C B is the transformation from the 
thruster reference frame to satellite reference frame, C B is the transformation from the satellite reference 
frame to ECI reference frame and [i is the Earth gravitational constant. In this paper, C B is chosen to be 
an identity matrix while it can be any other values. 

The attitude dynamics of the satellite is modelled for the nominal operation of the satellite where the 
three attitude angles Q = [</> 6 ip] T are defined as the angles between the nadir reference frame (i.e., the 
system of coordinates that maintain their orientation relative to the Earth as the satellite moves on its orbit) 
and satellite reference frame. The relationship between these two reference frames is 

W ?B = W NB + (16) 

where wf B is the angular velocity of satellite reference frame with respect to ECI reference frame expressed in 
the satellite reference frame, w^ B is the angular velocity of satellite reference frame with respect to the nadir 
reference frame expressed in the satellite reference frame, w^ N is the angular velocity of nadir reference frame 
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with respect to ECI reference frame expressed in the nadir reference frame and C ® is the transformation 
from the nadir reference frame to satellite reference frame. Assume that the attitude angles are small which 
is the case for an Earth pointing satellite. Then, w^ B can be approximated by 0 to the first order and from 
(16), 

0 = —C§w^ N + wf B (17) 


This can be used as the equation of motion associated with the attitude dynamics of the satellite given that 
wf B can be measured by the gyros. 

Therefore, by combining (15) and (17), the equation of motion of the satellite is 


d_ 

dt 


r 


r 

r 

= 

-^r + C^Cga 

© . 


. -C*w? N +wf B _ 


(18) 


where C B is the transformation from the satellite reference frame to nadir reference frame and C jy is the 
transformation from the nadir reference frame to ECI reference frame. Note that wf N and C'(- are functions 
of r and C B is a function of 0. Also note that a is the control input generated by the thrusters and iuf B is 
the measured input obtained by using the gyros. 


B. Measurement Model 

In this section, the measurement models of the star sensor, horizon sensor, sun sensor, GPS and inertial 
measurement unit are given. 


Star sensor Star sensor measures the azimuth < j> s tar and elevation A star of a star in the star sensor 
reference frame and provides attitude information of the satellite by comparing the measured star coordinate 
to the a priori known star coordinate in the star catalog. 20 The measured unit vector in the direction of the 
star in the star sensor reference frame is 


s = 


sin (f) s t ar COS A star 
COS (frstar COS A star 
Sill A star 


which can be related to the unit vector s in the direction of the star in ECI reference frame from the star 
catalog by 

8 = C s B tar C%C?s (19) 

where C B ar is the transformation from the satellite reference frame to star sensor reference frame. In this 
paper, C B ar is chosen to be an identity matrix and the star is chosen to be 15 degrees of right ascension and 
10 degrees of declination in ECI reference frame while they can be any other values. 


Horizon sensor Horizon sensor measures the angle between the satellite to Earth’s two horizons by 
scanning the Earth with the scanner’s line of sight sweeping from one horizon to the other and multiplying 
the time it takes to scan by the angular rate of the scanner rotation. 21 Horizon sensor usually operates 
with two synchronized scanners to scan upper and lower sides of the Earth that have equal distance to the 
equator of the Earth. Attitude information of the satellite can be obtained as follows. When there is no 
roll rotation, the times that it takes for the two scanners to scan the Earth are equal. When there is roll 
rotation, one scanner will take longer to scan the Earth than the other scanner. The difference of the time 
At can be related to the roll angle <f> by 


v At = 


4<p 


2 r 2 + R 2 e 
2 r 2 - R% 


(20) 


where v is the angular rate of the scanner rotation and R e is the radius of the Earth. When there is no pitch 
rotation, the satellite x-axis is perpendicular to the plane containing the two lines of sight at the middle of 
the scanning process. When there is pitch rotation, the angle between the satellite x-axis and the plane 
containing the two lines of sight at the middle of the scanning process is 90 degrees minus the pitch angle 6 , 


i.e., 



(21) 
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Sun SENSOR Sun sensor measures the azimuth 4> sun and elevation \ sun of sun in the sun sensor reference 
frame and provides attitude information of the satellite by comparing the measured sun coordinate to the 
actual sun coordinate. 20 The measured unit vector in the direction of sun in the sun sensor reference frame 
is 

I" COS (j) SU n sin A sun 1 


s = 


sin 4> sun sin A s 


COS A sun 


which can be related to the unit vector s in the direction of sun in ECI reference frame by 


£ /~isun s~iB s~iN 0 

h — Ug O * 


(22) 


where C s B n is the transformation from the satellite reference frame to sun sensor reference frame. In this 
paper, C^ n is chosen to be an identity matrix while it can be any other values. 


Other SENSORS For the GPS, the measurement is the position of the satellite in ECI reference frame, i.e., 


Vgps — r (23) 

For the inertial measurement unit, there are three gyros and three accelerometers. The three gyros measure 
the angular velocity of satellite reference frame with respect to ECI reference frame, i.e., 

Vgyro = w 1 b (24) 

The three accelerometers only measure the acceleration due to the three thrusters because they cannot mea- 
sure the gravitational force. By aligning the three accelerometers along the three thrusters, the measurement 
equation is 

y a = a (25) 


C. Linearized Satellite Model 

Because the fault detection filter design requires a linear model, the equation of motion (18) and measurement 
model are linearized around a Keplerian reference orbit with zero attitude angles. The Keplerian reference 
orbit is chosen with semimajor axis of 6998.12 km, eccentricity of 0.01, inclination of 5 deg, right ascension 
of the ascending node of 20 deg and argument of perigee of 15 deg. On the reference orbit, the nominal a is 
zero and the nominal wf B is equal to wf N . Let the state x and input u be 



r 


r -| 

X = 

r 

, u = 

a 


_ e _ 


Ugyro 


Denote the nominal state as xq and nominal input as uq. 
The linearized equation of motion of the satellite is 


Sr 



0 

I 


0 


Sr 


0 

0 ’ 


5a 

5ygyro 

Sr 

= 

dr 

dr 

Uo.MO 

0 


0 


Sr 

+ 

ril riN | 

\xo,uo 

0 


. SQ . 


de 

_ dr 

Uo>«o 

ae | 

Qrp \Xq,Uq 

00 

00 

Uo>«o 


. Se . 


0 

I 



(26) 


where some elements can be obtained analytically while other elements have to be obtained numerically 
using the central difference method because they are too complex. The linearized measurement model can 
be obtained similarly by linearizing the star sensor measurement y s tar — [ 4>star A stor ] , horizon sensor 
measurement yh = [v A t jh], sun sensor measurement y sun = [4>sun A stm ] and GPS measurement y gps 
analytically or numerically, i.e., 


6y = 


dystar 

Syh 

dysun 

dy gps 


= C 


Sr 

Sr 

56 


( 27 ) 
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IV. Periodic Fault Detection Filter Design and Evaluation 


In this section, periodic fault detection filters are designed and evaluated for the three thrusters, three 
gyros, star sensor, horizon sensor and sun sensor. Since there is redundancy in the number of GPS satellites, 
it is assumed that there is a separate health monitoring system built into the GPS receiver. 22 Therefore, it 
is assumed here that the GPS is always healthy. Since the three accelerometers measure only the control 
inputs but not the states, three algebraic parity equations can be constructed as 

r = y a ~ a 

For each residual, it is zero when there is no fault and becomes nonzero when the thruster fault or accelerom- 
eter fault occurs. Therefore, each parity equation can detect one thruster fault and one accelerometer fault, 
but cannot identify the fault that occurs. However, because the thruster fault can be detected and identified 
by using the periodic fault detection filter, the accelerometer fault can be detected and identified by using the 
parity equation assuming that these two faults do not occur at the same time. In Section A, periodic fault 
detection filters are designed and more algebraic parity equations are constructed. In Section B, periodic 
fault detection filters and parity equations are evaluated. 

A. Periodic Fault Detection Filter Design 

In order to design the periodic fault detection filter, the fault direction of each fault is obtained first. The 
fault directions for the three thruster faults are the first three columns of the B matrix. The fault directions 
for the three gyro faults are the last three columns of the B matrix. The fault directions for the star sensor, 
horizon sensor and sun sensor faults cannot be formed because the C matrix in (27) loses rank and there 
exists no f s in (3). Therefore, periodic fault detection filters are designed for the three thrusters and three 
gyros using (26) and (27) while periodic fault detection filters are designed for the star sensor, horizon sensor 
and sun sensor using (26) and a reduced-order measurement model that includes the GPS measurement and 
only one of the two components from the measurements of each of the star sensor, horizon sensor and sun 
sensor. 

For the three thruster and three gyro faults, six periodic single-fault detection filters are designed. For 
each single-fault detection filter, the target fault is one of the six faults and the nuisance fault is the other 
five faults. The filter gain of the periodic single-fault detection filter is given by (9) using the solution to 
(6). A process noise is assumed in (26) to model the uncertainty between the control commands and the 
control inputs generated by the thrusters and the sensor noise in the gyros. Two different sets of design 
weightings are chosen when designing the single-fault detection filters for the thruster faults. For the single- 
fault detection filters that detect the first and third thruster faults, the periodic Riccati matrix P,; in (6) is 
obtained with design weightings chosen as follows: Qi = 0.999, 7 = 10 — 2 , Qi = I, Q w is a diagonal matrix 
with 1, 1, 1, 5, 5, and 5 on the diagonal line and V is a diagonal matrix with 2.3373 x 10" 9 , 2.3373 x lO" 9 , 
7.2250 x 10 -5 , 7.6154 x 10~ 7 , 5.2716 x 10 -9 , and 5.2716 x 10 -9 on the diagonal line. For the single-fault 
detection filter that detects the second thruster fault, a larger Qi = 100/ is used to block the nuisance 
faults better. For the single-fault detection filters that detect the three gyro faults, the same set of design 
weightings is chosen. The periodic Riccati matrix Pi in (6) is obtained with design weightings chosen as 
follows: Qi = 0.99, 7 = 10 -2 , Qi = /, Q w = I and V = I. Note that V is treated as design parameter and 
chosen to be / instead of using the true sensor covariance when designing the single-fault detection filters 
for the three gyro faults because the performance of the single-fault detection filters is better when V = I. 

For the star sensor, horizon sensor and sun sensor faults, three periodic single-fault detection filters are 
designed similarly. For designing the single-fault detection filter that detects the star sensor fault, a reduced- 
order measurement model that includes y gps , Xstar, 7 h and <f> SU n is used. For designing the two single-fault 
detection filters that detect the horizon sensor and sun sensor faults, a reduced-order measurement model 
that includes y gps , 4>star, 7 h and X sun is used. These two reduced-order measurement models are used 
because they provide better rank condition than other six possible reduced-order measurement models over 
the entire period. 

The fault directions for the star sensor, horizon sensor and sun sensor faults in (3) are obtained by using 
fs = C T {CC T r'E s 

f s = C T (CC T )~ 1 E S - C T {CC T )-\CC T + CC t ){CC t )~ 1 E s 
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where C is the corresponding reduced-order measurement matrix. Note that / S ’s, the fault directions rep- 
resenting the fault magnitudes, are very small and about three orders smaller than / S ’s, the fault directions 
representing the fault rates. This indicates that if the faults are low-frequency, it will be difficult to detect 
the faults. For example, if the target fault is a bias, the projected residual will become small after a while 
even though the target fault still exists. Physically, this is because the attitude angles are weakly coupled 
with the rest of the dynamics and the star sensor, horizon sensor and sun sensor measure functions of mainly 
the attitude angles. On the other hand, since /„’ s are very small, the fault detection filter may be designed 
using only / S ’s. The advantage is that the fault detection filter can detect and identify more faults in one 
filter because the dimension of the fault is reduced. The disadvantage is that the projected residual may be 
a little more sensitive to the nuisance fault because the fault magnitude direction of the nuisance fault is not 
blocked. However, the fault magnitude direction of the nuisance fault may be included in the process noise 
to improve the blocking of the nuisance fault. 

When designing these three periodic single-fault detection filters, the target fault is one of the star sensor, 
horizon sensor and sun sensor faults and the nuisance fault is the other two sensor faults. A collection of the 
design weightings are used when tuning each of the single-fault detection filters. Then, the design weightings 
that offer the best performance are chosen. For the single-fault detection filter that detects the star sensor 
fault, the periodic Riccati matrix P,; in (6) is obtained with design weightings chosen as follows: Qi = 0.1, 
7 = 10 2 , Qi = /, V = I and Q w is a diagonal matrix with 1, 1, 1, 1, 1 and 10 4 on the diagonal line. 
Then, the filter gain of the periodic single-fault detection filter is obtained by (9). Similarly, the single-fault 
detection filter that detects the horizon sensor fault is designed with Qi = 0.85, 7 = 10 -2 , Qi = I, V = I and 
Q w = I. Finally, the single-fault detection filter that detects the sun sensor fault is designed with Q, = 0.1, 
7 = 10~ 2 , Qi is a diagonal matrix with 1 and 100 on the diagonal line, V = / and Q w is a diagonal matrix 
with 10 4 , 10 4 , 10 4 , 1, 10 2 and 10 4 on the diagonal line. 

Since the C matrix in (27) loses rank, there exists algebraic redundant relationship among the star sensor, 
horizon sensor, sun sensor and GPS that can be used to construct algebraic parity equations. First, there are 
three singular values of the C matrix that are very small compared to other singular values over the entire 
period. This indicates that there exist three parity equations for these sensors. In order to construct parity 
equations that are not sensitive to certain faults, three reduced-order measurement models, 5y r = C r Sx, are 
formed without using one of the star sensor, horizon sensor and sun sensor, respectively. For each of thee 
reduced-order measurement model, 5y r are the measurements from the GPS and two of the three sensors 
and 5x are the states defined in (27). For each of the three reduced-order measurement matrices C r , there 
is one singular value that is very small which indicates that there exists one parity equation which can be 
constructed as 

r = z T Sy r 

where z is the left singular vector of C r associated with the smallest singular value and z T C r = 0. Since 
each of these three parity equations is not sensitive to one of the star sensor, horizon sensor and sun sensor 
faults, these three parity equations can be used to detect and identify the star sensor, horizon sensor and 
sun sensor faults given that only one fault occurs at a time. 

B. Periodic Fault Detection Filter Evaluation 

These periodic fault detection filters and algebraic parity equations are implemented in simulation with 
measurements updated at 10 Hz. The parameters of the periodic fault detection filters and parity equations, 
and nominal measurement on the reference orbit are stored as a table with time of the reference orbit as the 
index over one period. In order to calculate the residuals of the fault detection filters and parity equations, 
the appropriate parameters and nominal measurement are obtained by using the point on the reference orbit 
that is closest to the position of the satellite as measured by GPS. Please refer to Ref. 14 for more details. 

The performance of the periodic fault detection filters and parity equations is evaluated when the satellite 
travels on another Keplerian orbit with the same orbital elements of the reference orbit except that the 
semimajor axis is 20 km larger. Realistic sensor noise is imposed on the sensors as white Gaussian noise 
with standard deviation of 4.8346 x 10 -5 rad for the star sensor, 0.0085 and 8.7266 x 10 -4 rad for the horizon 
sensor, 7.2606 x 10 -5 rad for the sun sensor, 1 m for the GPS and 0.0005 rad/s for the gyros. To evaluate 
the performance of the periodic fault detection filters and parity equations, a bias fault of 1 m/s 2 is imposed 
on each of the three thrusters separately at the 2000th second. Also, a bias fault of ten times larger than the 
standard deviation of the sensor noise is imposed on the star sensor, horizon sensor, sun sensor and three 
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gyros separately at the 2000th second. 

First, the performance of the six periodic single-fault detection filters for the three thrusters and three 
gyros is shown. Figure 1 shows the time response of the norm of the six projected residuals over one period 
when there is no fault. Figures 2 to 4 show the projected residuals when the first, second and third thruster 
fault occurs, respectively. Figures 5 to 7 show the projected residuals when the first, second and third gyro 
fault occurs, respectively. These figures show that when a thruster or gyro fault occurs, only the projected 
residual associated with the faulty thruster or gyro becomes large while the projected residuals associated 
with the healthy thrusters and gyros remain small. Therefore, the six periodic single-fault detection filters 
can detect and identify the three thruster faults and three gyro faults very well. Note that these projected 
residuals are scaled such that each projected residual becomes one when its associated target fault occurs 
and when the satellite travels on the reference orbit without any disturbances and uncertainty. 
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Figure 1. No fault occurs. 


Then, the performance of the three periodic single-fault detection filters and three parity equations for 
the star sensor, horizon sensor and sun sensor is shown. Figure 8 shows the time response of the norm of 
the three projected residuals of the periodic single-fault detection filters and the absolute value of the three 
residuals of the parity equations over one period when there is no fault. Figures 9 to 11 show the residuals 
when the star sensor, horizon sensor and sun sensor fault occurs, respectively. These figures show that 
when one of these three sensor faults occurs, only the projected residual of the periodic single-fault detection 
filter associated with the faulty sensor becomes large while the projected residuals of the periodic single-fault 
detection filters associated with the healthy sensors remain small. However, the projected residual associated 
with the faulty sensor becomes small after a while even though the fault still exists because the fault direction 
representing the fault magnitude is very small as discussed earlier. Note that the duration of the projected 
residual staying large is long enough for detecting and identifying the fault by using the fault detection filter 
along with the residual processor which is discussed in Section V. 

V. Residual Processor Design and Evaluation 

As discussed in the previous sections, there are nine periodic fault detection filters and three algebraic 
parity equations designed for residual generation. There are six elements from each filter residual and one 
form each parity equation. Therefore, the number of the input elements to the residual processor is 57 


10 of 22 


American Institute of Aeronautics and Astronautics 


First thruster residual 


Second thruster residual 


Third thruster residual 



Time (sec) 
First gyro residual 



1.5 


0.5 





Time (sec) 
Second gyro residual 


Time (sec) 
Third gyro residual 


Time (sec) 




Figure 2. First thruster fault occurs. 
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Figure 3. Second thruster fault occurs. 
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Figure 4. Third thruster fault occurs. 
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Figure 5. First gyro fault occurs. 
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Figure 6. Second gyro fault occurs. 
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Figure 7. Third gyro fault occurs. 
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Figure 8. No fault occurs. 
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Figure 9. Star sensor fault occurs. 
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Figure 10. Horizon sensor fault occurs. 
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Figure 11. Sun sensor fault occurs. 
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in total. This leads to intensive computations. To reduce the complexity of the problem and shorten the 
processing time, the norm of the residual of each filter and parity equation is used instead. Then, the input 
vector of the residual processor has only 12 elements. Another advantage of using just the norm is that only 
one hypothesis is enough when defining a failure hypothesis. Without the norm operating on the residual, 
two hypotheses are needed to define one fault. One hypothesis is needed to define the positive fault scenario 
and the another hypothesis is associated to the negative fault scenario. 

Including null hypothesis and nine failure hypotheses, there are totally ten hypotheses are defined when 
designing the residual processor. In tradition, the no fault hypothesis is referred to as the null hypothesis. 
The nine failure hypotheses are first thruster fault hypothesis, second thruster fault hypothesis, third thruster 
fault hypothesis, first gyro fault hypothesis, second gyro fault hypothesis, third gyro fault hypothesis, star 
sensor fault hypothesis, horizon sensor fault hypothesis and sun sensor fault hypothesis. The output of the 
residual processor are the probabilities of the ten hypotheses conditioned on the history of the residuals. 

Now, we are ready to define the hypothesis conditioned density functions. For the problem at hand, 
it is assumed that the input vector generated at each time instant is Gaussian distributed. Following the 
formulation shown in (14), a Gaussian distributed random vector x £ lZ n where x ~ J\f(rrii,A x .) with the 
mean rrq £ TZ n and to* ~ Unif[bi , bi + 2 • m*] Vto* > 0 where bi £ 1Z n , the probability density function can 
be written as 


fi(x) = 




I er /{ V5 ’ Ax ' /2 ' ^ ~ bi ^ ~ er/{ 7f ' A ** 1/2 ’ 0 “ b i ~ 2 ’ m i)}] 


Note that Gaussian model may not be the most accurate assumption. Other probability density functions 
can be chosen to propagate the posterior probability. The analysis in the Shiryayev test still remains the 
same. 

The mean and covariance of the conditional probability density under each of the given hypothesis 
are determined from the simulation. For example, if the magnitude of all the filter residuals are smaller 
than 0.5 when no fault presents, the mean of the input residual vector can be assumed to be uniformly 
distributed between 0 and 0.5 when defining the null hypothesis. To define a particular failure hypothesis, 
the corresponding fault with a variety of magnitude are introduced in the simulation. Then, the corresponding 
magnitude of the residuals are recorded and the range of the mean can be determined accordingly. Finally, 
the covariances can be calculated from the simulated residuals. In our analysis, the priori probability of 
transition pi is assumed to be ICG 8 and the a priori probabilities 7r,; is chosen to be 0.001. Both of them are 
assumed to be same for all hypotheses. 

The multiple-hypothesis Shiryayev sequential probability test can now be applied on the residuals of 
the periodic fault detection filters and the parity equations to propagate the posterior probability of each 
hypothesis. Figure 12 shows the posterior probability conditioned on each of the nine failure hypothesis when 
there is no fault. It can be seen that the failure probabilities stay at zero in the entire period. Figure 13 
to 15 shows the hypothesis conditioned probability when the first, second and third thruster fault occurs, 
respectively. It is shown that our scheme correctly identified the three thruster faults. Figure 16 to 18 show 
the conditioned probability when the first, second and third gyro fault occurs, respectively. The processor 
has a little trouble in differentiating between first gyro fault and second gyro fault right after the fault 
is occurred. However, this should not cause a false declaration when the false alarm rate is chosen to be 
smaller than 0.1 which leads to a threshold probability bigger than 0.9. Therefore, the Shiryayev test can still 
correctly identifies these two failures. Finally, Figures 19 to 21show the posterior probability when the star, 
horizon and sun sensor fault occurs, respectively. The results show that the sensor faults can be detected 
and identified within fraction of a second. 


VI. Conclusion 

A complete health monitoring system based on analytical redundancy is developed for satellites on 
elliptical orbits which can be modelled as periodic systems. This health monitoring system consists of 
periodic fault detection filters, algebraic parity equations and residual processor. The actuators and sensors 
being monitored include three thrusters, three gyros, a star sensor, a horizon sensor, a sun sensor and 
three accelerometers. The performance of the health monitoring system is evaluated in the simulation in 
the presence of disturbances and uncertainty. For the star sensor, horizon sensor and sun sensor faults, 
the residual processor announces the fault in a fraction of second. For the three gyro faults, the residual 
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Figure 12. No fault occurs. 
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Figure 13. First thruster fault occurs. 
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Figure 14. Second thruster fault occurs. 
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Figure 15. Third thruster fault occurs. 
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Figure 16. 


First gyro fault occurs. 
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Figure 17. Second gyro fault occurs. 
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Figure 18. Third gyro fault occurs. 


First thruster hypothesis 


0.5 


Second thruster hypothesis Third thruster hypothesis 

1 


0.5 


0.5 


1999 2000 2001 2002 1999 2000 2001 2002 1999 2000 2001 2002 


First gyro hypothesis 


Second gyro hypothesis 


Third gyro hypothesis 


0.5 


0 



1999 2000 2001 2002 1999 2000 2001 2002 1999 2000 2001 2002 


Star sensor hypothesis 


Horizon sensor hypothesis 


Sun sensor hypothesis 



0.5 


0.5 


Time (sec) Time (sec) Time (sec) 

Figure 19. Star sensor fault occurs. 
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Figure 20. Horizon sensor fault occurs. 
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Figure 21. Sun sensor fault occurs. 
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processor needs up to about one second to announce the fault. For the three thruster faults, the residual 
processor needs about four seconds to announce the fault. Note that the satellite will only deviate from the 
desired trajectory by 8 to in four seconds given that the thruster fault is 1 m/s 2 . When the fault magnitude 
becomes larger, the time that the residual processor needs to announce the fault will decrease. In this paper, 
only single satellite is considered. When a cluster of satellites are considered, the time that the residual 
processor needs to detect the thruster fault is expected to decrease because the additional range sensor can 
be very accurate. 
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